This knitr document contains the code used to analyze and visualize data from the VIS camera system and water data. Knit-ing this document will generate an HTML document that contains both the embedded R code and the output.

Attach the required R packages

library(qvalue, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)
library(lubridate, warn.conflicts = FALSE)
library(MASS, warn.conflicts = FALSE)
library(car, warn.conflicts = FALSE)
## Warning: package 'car' was built under R version 3.1.2
library(nlme, warn.conflicts = FALSE)
library(mvtnorm, warn.conflicts = FALSE)
library(grid, warn.conflicts = FALSE)

Soil volume water content

Use linear regression to create a simple model for using water volume to predict soil volume water content.

Download data for soil wet/dry weight and volume water content measurements.

if (!file.exists('soil_weigth_vwc.txt')) {
  download.file('http://files.figshare.com/1939954/soil_weigth_vwc.txt',
                'soil_weigth_vwc.txt')
}

Read the soil volume water content data

vwc.data = read.table(file="soil_weigth_vwc.txt", sep="\t", header=TRUE)

Calculate the average soil dry weight

mean(vwc.data$weight_dry)
## [1] 72.92138

Create a linear model for water volume and volume water content. Water volumes >= 260 appear to have saturated the soil water carrying capacity, so remove them from the model.

vwc.lm = lm(vwc_wet ~ water_vol, data=vwc.data[vwc.data$water_vol < 260,])
summary(vwc.lm)
## 
## Call:
## lm(formula = vwc_wet ~ water_vol, data = vwc.data[vwc.data$water_vol < 
##     260, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5426 -1.6196  0.3976  1.4214  4.0771 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.032863   1.389638  -2.182   0.0391 *  
## water_vol    0.233197   0.008145  28.630   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.082 on 24 degrees of freedom
## Multiple R-squared:  0.9716, Adjusted R-squared:  0.9704 
## F-statistic: 819.7 on 1 and 24 DF,  p-value: < 2.2e-16

Plot the model

vwc.plot = ggplot(vwc.data[vwc.data$water_vol < 260,], aes(x=water_vol, y=vwc_wet)) +
                  geom_point(size=2) +
                  geom_smooth(method='lm', formula=y~x) +
                  scale_x_continuous("Water volume (ml)") +
                  scale_y_continuous("Volume water content (%)") +
                  theme_bw() +
                  theme(axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold"))
print(vwc.plot)

Predict volume water contents for the water treatment groups

treatment.water = data.frame(water_vol=c(217, 144.5, 72))
treatment.water$vwc = predict(object = vwc.lm, newdata=treatment.water)
print(treatment.water)
##   water_vol      vwc
## 1     217.0 47.57084
## 2     144.5 30.66407
## 3      72.0 13.75731

Convert VIS camera zoom units

LemnaTec VIS camera zoom units range from 1 to 6000, which correspond to 1 to 6X zoom.

zoom.lm = lm(zoom.camera ~ zoom, data=data.frame(zoom=c(1,6000),
                                                    zoom.camera=c(1,6)))
summary(zoom.lm)
## 
## Call:
## lm(formula = zoom.camera ~ zoom, data = data.frame(zoom = c(1, 
##     6000), zoom.camera = c(1, 6)))
## 
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9991665         NA      NA       NA
## zoom        0.0008335         NA      NA       NA
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 1 and 0 DF,  p-value: NA

VIS zoom correction

In this section we define models that are used to convert area and length between camera zoom levels to a common scale.

Download data for a reference object imaged at different zoom levels.

if (!file.exists('zoom_calibration_data.txt')) {
  download.file('http://files.figshare.com/1845355/zoom_calibration_data.txt',
                'zoom_calibration_data.txt')
}

Read zoom calibrartion data

z.data = read.table(file="zoom_calibration_data.txt", sep="\t", header=TRUE)

Calculate px per cm

z.data$px_cm = z.data$length_px / z.data$length_cm

Convert LemnaTec zoom units to camera zoom units

z.data$zoom.camera = predict(object = zoom.lm, newdata=z.data)

Zoom correction for area

Fit a variety of regression models to relative object area by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC)

Non-linear regression (exponential)

area.coef = coef(nls(log(rel_area) ~ log(a * exp(b * zoom.camera)),
                     z.data, start = c(a = 1, b = 0.01)))
area.coef = data.frame(a=area.coef[1], b=area.coef[2])
area.nls = nls(rel_area ~ a * exp(b * zoom.camera),
               data = z.data, start=c(a=area.coef$a, b=area.coef$b))
summary(area.nls)
## 
## Formula: rel_area ~ a * exp(b * zoom.camera)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a  0.37245    0.02159   17.25   <2e-16 ***
## b  0.86668    0.01504   57.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3306 on 32 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 1.605e-06

Non-linear regression (polynomial)

area.pol = lm(rel_area ~ zoom.camera + I(zoom.camera^2), z.data)
summary(area.pol)
## 
## Call:
## lm(formula = rel_area ~ zoom.camera + I(zoom.camera^2), data = z.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8494 -0.3820  0.1176  0.3173  0.7494 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.42793    0.48203   9.186 2.34e-10 ***
## zoom.camera      -4.68583    0.40917 -11.452 1.14e-12 ***
## I(zoom.camera^2)  1.64781    0.07786  21.163  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.431 on 31 degrees of freedom
## Multiple R-squared:  0.9907, Adjusted R-squared:   0.99 
## F-statistic:  1643 on 2 and 31 DF,  p-value: < 2.2e-16

AIC

AIC(area.nls, area.pol)
##          df      AIC
## area.nls  3 25.16270
## area.pol  4 44.12291

The exponential model minimizes AIC. Plot exponential model.

nls.plot = ggplot(z.data, aes(x=zoom.camera, y=rel_area)) +
                  geom_point(size=2.5) +
                  scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
                  scale_y_continuous(lim=c(0,16),
                                     "Reference object relative pixel area") +
                  stat_smooth(data=z.data, aes(x=zoom.camera, y=rel_area),
                              method="nls", se=FALSE, formula=y ~ a * exp(b * x),
                              start=c(a=area.coef$a, b=area.coef$b)) +
                  theme_bw() +
                  theme(axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold"))
nls.plot = nls.plot + labs(title='Figure 12A')
print(nls.plot)

Zoom correction for length

Fit a variety of regression models to px/cm by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC).

Non-linear regression (exponential)

len.coef = coef(nls(log(px_cm) ~ log(a * exp(b * zoom.camera)),
                    z.data[z.data$camera == 'VIS SV',], start = c(a = 1, b = 0.01)))
len.coef = data.frame(a=len.coef[1], b=len.coef[2])
len.nls = nls(px_cm ~ a * exp(b * zoom.camera),
              data = z.data[z.data$camera == 'VIS SV',],
              start=c(a=len.coef$a, b=len.coef$b))
summary(len.nls)
## 
## Formula: px_cm ~ a * exp(b * zoom.camera)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 8.370121   0.150430   55.64   <2e-16 ***
## b 0.445769   0.005504   80.98   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5106 on 14 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 2.15e-06

Length zoom correction using a 2 order polynomial (px/cm given a zoom setting). Length correction only works for side-view images right now because scale in top-view images are affected by the plant lifter position in addition to zoom.

len.poly = lm(px_cm ~ zoom.camera + I(zoom.camera^2),
              data=z.data[z.data$camera == 'VIS SV',])
summary(len.poly)
## 
## Call:
## lm(formula = px_cm ~ zoom.camera + I(zoom.camera^2), data = z.data[z.data$camera == 
##     "VIS SV", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29376 -0.11153  0.02836  0.14223  0.21043 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      14.89815    0.30540   48.78 4.13e-16 ***
## zoom.camera      -3.74980    0.27285  -13.74 4.04e-09 ***
## I(zoom.camera^2)  3.12322    0.05472   57.08  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1741 on 13 degrees of freedom
## Multiple R-squared:  0.9998, Adjusted R-squared:  0.9998 
## F-statistic: 3.402e+04 on 2 and 13 DF,  p-value: < 2.2e-16

AIC

AIC(len.nls, len.poly)
##          df       AIC
## len.nls   3 27.762784
## len.poly  4 -5.857059

The polynomial model minimizes AIC. Plot polynomial model.

poly.plot = ggplot(z.data[z.data$camera == 'VIS SV',], aes(x=zoom.camera, y=px_cm)) +
                   geom_point(size=2.5) +
                   scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
                   scale_y_continuous(lim=c(0,50),
                                      "Reference object length scale (px/cm)") +
                   stat_smooth(data=z.data[z.data$camera == 'VIS SV',],
                               aes(x=zoom.camera, y=px_cm), method="lm",
                               formula=y ~ x + I(x^2)) +
                   theme_bw() +
                   theme(axis.title.x=element_text(face="bold"),
                         axis.title.y=element_text(face="bold"))
poly.plot = poly.plot + labs(title='Figure 12B')
print(poly.plot)

Height modeling

In this section we measure how well PlantCV estimates plant height.

Download data manually measured plant height data set (n = 173).

if (!file.exists('height_model_data.txt')) {
  download.file('http://files.figshare.com/1845361/height_model_data.txt',
                'height_model_data.txt')
}

Read height data.

ht.data = read.table(file="height_model_data.txt",sep="\t",header=TRUE)

Convert LemnaTec zoom units to camera zoom units

ht.data$zoom.camera = predict(object = zoom.lm, newdata=ht.data)

Convert pixel traits to cm with zoom correction. px/cm calibration

px_cm = predict(object = len.poly, newdata=ht.data)
ht.data$height_above_bound_cm = ht.data$height_above_bound / px_cm

Height linear model

height.model = lm(manual_height~height_above_bound_cm, ht.data)
summary(height.model)
## 
## Call:
## lm(formula = manual_height ~ height_above_bound_cm, data = ht.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9229 -0.3034 -0.1263  0.3916  5.2015 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.303813   0.149776   2.028   0.0441 *  
## height_above_bound_cm 0.852158   0.002789 305.528   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.103 on 171 degrees of freedom
## Multiple R-squared:  0.9982, Adjusted R-squared:  0.9982 
## F-statistic: 9.335e+04 on 1 and 171 DF,  p-value: < 2.2e-16

Plot the height linear model.

hlm.plot = ggplot(ht.data, aes(x=height_above_bound_cm, y=manual_height)) +
                  geom_point(aes(color=factor(round(zoom.camera,1))),size=2.5) +
                  geom_smooth(method="lm", formula=y~x, color='black') +
                  scale_x_continuous(lim=c(0,110), breaks=c(0,25,50,75,100), "Estimated height (cm)") +
                  scale_y_continuous(lim=c(0,110), breaks=c(0,25,50,75,100), "Manually measured height (cm)") +
                  theme_bw() +
                  theme(legend.position=c(0.2,0.75),
                        axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold")) +
                  labs(color="Zoom")
hlm.plot = hlm.plot + labs(title='Figure 3A')
print(hlm.plot)

Biomass modeling

In this section we model fresh- and dry-weight above ground biomass using image measurements.

Download data manually measured plant biomass data set (n = 41).

if (!file.exists('biomass_model_data.txt')) {
  download.file('http://files.figshare.com/1845360/biomass_model_data.txt',
                'biomass_model_data.txt')
}

Read biomass data.

st.data = read.table(file='biomass_model_data.txt', sep="\t", header=TRUE,
                     stringsAsFactors=FALSE)

Create out-of-frame indicator variable.

st.data$outind = NA
st.data[st.data$out_of_frame == T,]$outind = 1
st.data[st.data$out_of_frame == F,]$outind = 0

Genotype indicator variable.

st.data$group = NA
st.data[st.data$genotype == 'A10',]$group = 0
st.data[st.data$genotype == 'B100',]$group = 1

Fresh-weight biomass

Full model. Includes side-view area, top-view area and height.

fw.full = lm(fresh_weight ~ sv_area * tv_area * height_above_bound, st.data)

Step-wise model selection with AIC.

fw.step = stepAIC(fw.full, direction="both")
## Start:  AIC=-50.4
## fresh_weight ~ sv_area * tv_area * height_above_bound
## 
##                                      Df Sum of Sq    RSS     AIC
## - sv_area:tv_area:height_above_bound  1  0.055428 8.1731 -52.122
## <none>                                            8.1177 -50.401
## 
## Step:  AIC=-52.12
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area + 
##     sv_area:height_above_bound + tv_area:height_above_bound
## 
##                                      Df Sum of Sq    RSS     AIC
## - tv_area:height_above_bound          1  0.024335 8.1974 -54.000
## - sv_area:tv_area                     1  0.061951 8.2350 -53.812
## - sv_area:height_above_bound          1  0.108786 8.2819 -53.580
## <none>                                            8.1731 -52.122
## + sv_area:tv_area:height_above_bound  1  0.055428 8.1177 -50.401
## 
## Step:  AIC=-54
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area + 
##     sv_area:height_above_bound
## 
##                              Df Sum of Sq    RSS     AIC
## - sv_area:tv_area             1   0.31153 8.5089 -54.471
## <none>                                    8.1974 -54.000
## - sv_area:height_above_bound  1   0.68121 8.8786 -52.727
## + tv_area:height_above_bound  1   0.02433 8.1731 -52.122
## 
## Step:  AIC=-54.47
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound
## 
##                              Df Sum of Sq    RSS     AIC
## <none>                                    8.5089 -54.471
## - sv_area:height_above_bound  1   0.42578 8.9347 -54.469
## + sv_area:tv_area             1   0.31153 8.1974 -54.000
## - tv_area                     1   0.54733 9.0563 -53.915
## + tv_area:height_above_bound  1   0.27391 8.2350 -53.812
summary(fw.step)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound + 
##     sv_area:height_above_bound, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.19711 -0.10790 -0.00974  0.08153  1.55179 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.885e-01  2.513e-01   0.750   0.4581    
## sv_area                     4.673e-05  4.812e-06   9.710 1.36e-11 ***
## tv_area                    -1.045e-05  6.864e-06  -1.522   0.1368    
## height_above_bound         -3.280e-02  1.451e-02  -2.260   0.0299 *  
## sv_area:height_above_bound  1.166e-07  8.688e-08   1.342   0.1879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4862 on 36 degrees of freedom
## Multiple R-squared:  0.9835, Adjusted R-squared:  0.9816 
## F-statistic: 535.9 on 4 and 36 DF,  p-value: < 2.2e-16

AIC model

fw.aic = lm(fresh_weight ~ sv_area + tv_area + height_above_bound +
              sv_area*height_above_bound, st.data)
summary(fw.aic)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound + 
##     sv_area * height_above_bound, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.19711 -0.10790 -0.00974  0.08153  1.55179 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.885e-01  2.513e-01   0.750   0.4581    
## sv_area                     4.673e-05  4.812e-06   9.710 1.36e-11 ***
## tv_area                    -1.045e-05  6.864e-06  -1.522   0.1368    
## height_above_bound         -3.280e-02  1.451e-02  -2.260   0.0299 *  
## sv_area:height_above_bound  1.166e-07  8.688e-08   1.342   0.1879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4862 on 36 degrees of freedom
## Multiple R-squared:  0.9835, Adjusted R-squared:  0.9816 
## F-statistic: 535.9 on 4 and 36 DF,  p-value: < 2.2e-16

The AIC model contains tv_area and height which does not have a significant coefficient, test dropping.

fw.red = lm(fresh_weight ~ sv_area, st.data)
summary(fw.red)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17949 -0.25179  0.02035  0.23804  1.28235 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.044e-01  1.060e-01  -2.872  0.00657 ** 
## sv_area      3.935e-05  8.757e-07  44.931  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5003 on 39 degrees of freedom
## Multiple R-squared:  0.981,  Adjusted R-squared:  0.9806 
## F-statistic:  2019 on 1 and 39 DF,  p-value: < 2.2e-16

Goodness of fit.

anova(fw.aic, fw.red)
## Analysis of Variance Table
## 
## Model 1: fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area * 
##     height_above_bound
## Model 2: fresh_weight ~ sv_area
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     36 8.5089                           
## 2     39 9.7630 -3   -1.2541 1.7686 0.1706

SV area model.

sv.model = lm(fresh_weight ~ sv_area, st.data)
summary(sv.model)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17949 -0.25179  0.02035  0.23804  1.28235 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.044e-01  1.060e-01  -2.872  0.00657 ** 
## sv_area      3.935e-05  8.757e-07  44.931  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5003 on 39 degrees of freedom
## Multiple R-squared:  0.981,  Adjusted R-squared:  0.9806 
## F-statistic:  2019 on 1 and 39 DF,  p-value: < 2.2e-16

Plot SV model

sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=fresh_weight)) +
                       geom_smooth(method="lm", color="black", formula = y ~ x) +
                       geom_point(size=2.5) +
                       scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                       scale_y_continuous("Fresh-weight biomass (g)") +
                       theme_bw() +
                       theme(axis.title.x=element_text(face="bold"),
                             axis.title.y=element_text(face="bold"))
sv.model.plot = sv.model.plot + labs(title='Figure 4A')
print(sv.model.plot)

SV area model with out-of-frame.

sv.ind.model = lm(fresh_weight ~ sv_area + outind, st.data)
anova(sv.model, sv.ind.model)
## Analysis of Variance Table
## 
## Model 1: fresh_weight ~ sv_area
## Model 2: fresh_weight ~ sv_area + outind
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)   
## 1     39 9.7630                              
## 2     38 7.5498  1    2.2132 11.139 0.0019 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(sv.ind.model)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area + outind, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75671 -0.21287 -0.00372  0.22304  1.17010 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.499e-01  9.583e-02  -2.607   0.0130 *  
## sv_area      4.028e-05  8.288e-07  48.599   <2e-16 ***
## outind      -5.963e-01  1.787e-01  -3.338   0.0019 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4457 on 38 degrees of freedom
## Multiple R-squared:  0.9853, Adjusted R-squared:  0.9846 
## F-statistic:  1277 on 2 and 38 DF,  p-value: < 2.2e-16

Plot SV model with out-of-frame.

sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight, 
                                        group=out_of_frame, color=out_of_frame)) +
                           geom_point(size=2.5) +
                           geom_smooth(method="lm", color="black") +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Fresh-weight biomass (g)") +
                           theme_bw() +
                           theme(legend.position=c(0.2,0.8),
                                 axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold")) +
                           labs(color="Out-of-frame")
print(sv.model.out.plot)

SV area model with genotype

sv.gt.model = lm(fresh_weight ~ sv_area + group, st.data)
summary(sv.gt.model)
## 
## Call:
## lm(formula = fresh_weight ~ sv_area + group, data = st.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.32109 -0.29602  0.03093  0.21960  1.15107 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.510e-01  1.342e-01  -3.362  0.00178 ** 
## sv_area      3.956e-05  8.637e-07  45.802  < 2e-16 ***
## group        2.647e-01  1.542e-01   1.717  0.09419 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4883 on 38 degrees of freedom
## Multiple R-squared:  0.9824, Adjusted R-squared:  0.9815 
## F-statistic:  1061 on 2 and 38 DF,  p-value: < 2.2e-16

Plot SV model with genotype.

sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight, 
                                        group=genotype, color=genotype)) +
                           geom_point(size=2.5) +
                           geom_smooth(method="lm", color="black") +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Fresh-weight biomass (g)") +
                           theme_bw() +
                           theme(legend.position=c(0.2,0.8),
                                 axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold")) +
                           labs(color="Genotype")
print(sv.model.gen.plot)

Dry-weight biomass

SV area model

dry.sv.model = lm(dry_weight ~ sv_area, st.data)
summary(dry.sv.model)
## 
## Call:
## lm(formula = dry_weight ~ sv_area, data = st.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.183938 -0.047916  0.009925  0.039315  0.195700 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.295e-02  1.479e-02   -3.58 0.000938 ***
## sv_area      4.634e-06  1.222e-07   37.93  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06981 on 39 degrees of freedom
## Multiple R-squared:  0.9736, Adjusted R-squared:  0.9729 
## F-statistic:  1438 on 1 and 39 DF,  p-value: < 2.2e-16

Plot dry-weight side-view model

dry.sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=dry_weight)) +
                           geom_smooth(method="lm", color="black", formula = y ~ x) +
                           geom_point(size=2.5) +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Dry-weight biomass (g)") +
                           theme_bw() +
                           theme(axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold"))
dry.sv.model.plot = dry.sv.model.plot + labs(title='Supplemental Figure S2')
print(dry.sv.model.plot)

SV area model with out-of-frame.

dry.sv.ind.model = lm(dry_weight ~ sv_area + outind, st.data)
summary(dry.sv.ind.model)
## 
## Call:
## lm(formula = dry_weight ~ sv_area + outind, data = st.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.133525 -0.027886  0.006673  0.024790  0.131393 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.393e-02  1.255e-02  -3.501 0.001202 ** 
## sv_area      4.789e-06  1.085e-07  44.124  < 2e-16 ***
## outind      -9.869e-02  2.340e-02  -4.219 0.000147 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05836 on 38 degrees of freedom
## Multiple R-squared:  0.982,  Adjusted R-squared:  0.9811 
## F-statistic:  1038 on 2 and 38 DF,  p-value: < 2.2e-16
dry.sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight, 
                                        group=out_of_frame, color=out_of_frame)) +
                           geom_point(size=2.5) +
                           geom_smooth(method="lm", color="black") +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Dry-weight biomass (g)") +
                           theme_bw() +
                           theme(legend.position=c(0.2,0.8),
                                 axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold")) +
                           labs(color="Out-of-frame")
print(dry.sv.model.out.plot)

SV area model with genotypes

dry.sv.gt.model = lm(dry_weight ~ sv_area + group, st.data)
summary(dry.sv.gt.model)
## 
## Call:
## lm(formula = dry_weight ~ sv_area + group, data = st.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.183197 -0.049205  0.008654  0.040345  0.196263 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.413e-02  1.943e-02  -2.786  0.00828 ** 
## sv_area      4.636e-06  1.251e-07  37.061  < 2e-16 ***
## group        2.130e-03  2.233e-02   0.095  0.92451    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07071 on 38 degrees of freedom
## Multiple R-squared:  0.9736, Adjusted R-squared:  0.9722 
## F-statistic:   701 on 2 and 38 DF,  p-value: < 2.2e-16
dry.sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight, 
                                        group=genotype, color=genotype)) +
                           geom_point(size=2.5) +
                           geom_smooth(method="lm", color="black") +
                           scale_x_continuous("Shoot and leaf area (x10^5 px)") +
                           scale_y_continuous("Dry-weight biomass (g)") +
                           theme_bw() +
                           theme(legend.position=c(0.2,0.8),
                                 axis.title.x=element_text(face="bold"),
                                 axis.title.y=element_text(face="bold")) +
                           labs(color="Genotype")
print(dry.sv.model.gen.plot)

Analyze VIS data

In this section we analyze geometric traits from the complete VIS data set.

Download the complete VIS data set. There were a total of 6,399 snapshots with VIS image data, but the download only includes the 6,207 snapshots that were successfully processed by PlantCV. Failed snapshots generally are those that lack a plant (empty pot controls, dead plants, etc.)

if (!file.exists('vis_snapshots.txt')) {
  download.file('http://files.figshare.com/1845362/vis_snapshots.txt',
                'vis_snapshots.txt')
}

Read data and format for analysis

vis.data = read.table(file="vis_snapshots.txt", sep=",", header=TRUE)

Planting date

planting_date = as.POSIXct("2013-11-26")

Add water treatment column coded in barcodes.

vis.data$treatment <- NA
vis.data$treatment[grep("AA", vis.data$plant_id)] <- 100
vis.data$treatment[grep("AB", vis.data$plant_id)] <- 0
vis.data$treatment[grep("AC", vis.data$plant_id)] <- 16
vis.data$treatment[grep("AD", vis.data$plant_id)] <- 33
vis.data$treatment[grep("AE", vis.data$plant_id)] <- 66

Add plant genotype column coded in barcodes.

vis.data$genotype <- NA
vis.data$genotype[grep("p1", vis.data$plant_id)] <- 'A10'
vis.data$genotype[grep("p2", vis.data$plant_id)] <- 'B100'
vis.data$genotype[grep("r1", vis.data$plant_id)] <- 'R20'
vis.data$genotype[grep("r2", vis.data$plant_id)] <- 'R70'
vis.data$genotype[grep("r3", vis.data$plant_id)] <- 'R98'
vis.data$genotype[grep("r4", vis.data$plant_id)] <- 'R102'
vis.data$genotype[grep("r5", vis.data$plant_id)] <- 'R128'
vis.data$genotype[grep("r6", vis.data$plant_id)] <- 'R133'
vis.data$genotype[grep("r7", vis.data$plant_id)] <- 'R161'
vis.data$genotype[grep("r8", vis.data$plant_id)] <- 'R187'

Add genotype x treatment group column

vis.data$group = paste(vis.data$genotype,'-',vis.data$treatment,sep='')

Remove plants that were sampled for biomass measurements.

vis.data = vis.data[!vis.data$plant_id %in% st.data$barcode,]

Add calendar-time data column using the Unix-time data

vis.data$date = as.POSIXct(vis.data$datetime, origin = "1970-01-01")

Calculate days after planting from planting data

vis.data$dap = as.numeric(vis.data$date - planting_date)

Predict fresh- and dry-weight biomass from linear models

vis.data$fw_biomass = predict.lm(object = sv.model, newdata=vis.data)
vis.data$dw_biomass = predict.lm(object = dry.sv.model, newdata=vis.data)

Plant height analysis

Plot height for S. viridis and S. italica water treatments 100% and 33% full-capacity.

height.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
                               vis.data$genotype == 'B100') &
                              (vis.data$treatment == 100 |
                               vis.data$treatment == 33),],
                     aes(x=dap, y=height_above_bound, color=factor(group))) +
                     geom_point(size=2.5) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(0,120),
                                        breaks=c(0,20,40,60,80,100,120),
                                        name="Estimated height (cm)") +
                     theme_bw() +
                     theme(legend.position=c(0.25,0.75),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Genotype-treatment")
height.plot = height.plot + labs(title='Figure 3B')
print(height.plot)

Plot height for all genotypes.

height.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
                                    vis.data$treatment == 33,],
                           aes(x=dap, y=height_above_bound, color=factor(treatment))) +
                     geom_point(size=1) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(0,120),
                                        breaks=c(0,20,40,60,80,100,120),
                                        name="Estimated height (cm)") +
                     theme_bw() +
                     theme(legend.position=c(0.8,0.1),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Treatment") +
                     facet_wrap(~ genotype, ncol = 3)
height.facet.plot = height.facet.plot + labs(title='Supplemental Figure S1')
print(height.facet.plot)

Genotype differences in height

For each day, calculate the mean and 95% confidence intervals for each genotype in the control water group.

height_per_day = function(vis.data) {
  dap = c()
  genotype = c()
  int.low = c()
  int.up = c()
  est = c()
  
  # Loop through each day and genotype
  for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
    if(d != 24) {
      for(g in levels(factor(vis.data$genotype))) {
        h.data = vis.data[vis.data$genotype == g & vis.data$treatment == 100 &
                          as.integer(vis.data$dap) == d,]$height_above_bound
        if(sd(h.data) != 0) {
          dap = c(dap,d)
          genotype = c(genotype,g)
          test = t.test(h.data)
          int.low = c(int.low,test$conf.int[1])
          int.up = c(int.up,test$conf.int[2])
          est = c(est,test$estimate)
        }
      }
    }
  }
  #drought_resp = drought_resp[-14,]
  results = data.frame(dap=as.numeric(dap),
                       genotype=genotype,
                       conf.int.low=int.low,
                       conf.int.up=int.up,
                       mean=est)  

  return(results)
}

Height per day 95% CI. Write results to file.

height.perday.results = height_per_day(vis.data)
write.csv(height.perday.results,file="height.perday.100water.csv")
print(height.perday.results)
##     dap genotype conf.int.low conf.int.up       mean
## 1    11      A10    7.8384656    9.927418   8.882942
## 2    11     B100    5.3761414    6.677522   6.026831
## 3    11     R102    2.6078989    6.747923   4.677911
## 4    11     R128    3.4255116    6.098742   4.762127
## 5    11     R133    5.8106341   11.190037   8.500336
## 6    11     R161    5.9550014    6.742792   6.348897
## 7    11     R187    4.6001185    8.238236   6.419177
## 8    11      R20    8.7685119   12.487009  10.627761
## 9    11      R70    6.9589765    8.860251   7.909614
## 10   11      R98   11.3958077   13.499331  12.447569
## 11   12      A10   10.5706419   14.040344  12.305493
## 12   12     B100    7.1779957    9.543168   8.360582
## 13   12     R102    4.1386602   10.394107   7.266383
## 14   12     R128    4.7644938    9.358302   7.061398
## 15   12     R133    0.3769178   32.979175  16.678046
## 16   12     R161    5.5667113   10.395653   7.981182
## 17   12     R187    8.9117000   14.337900  11.624800
## 18   12      R20  -24.8325684   57.616115  16.391774
## 19   12      R70    6.2768196   11.896208   9.086514
## 20   12      R98   14.4614239   19.097635  16.779529
## 21   13      A10   14.2175454   16.845714  15.531630
## 22   13     B100   10.9868413   12.771688  11.879265
## 23   13     R102    8.1006537   12.638227  10.369440
## 24   13     R128   10.0067593   13.083564  11.545162
## 25   13     R133   15.6105089   20.434429  18.022469
## 26   13     R161   11.5326626   14.465157  12.998910
## 27   13     R187   11.0472672   15.088387  13.067827
## 28   13      R20   14.8393242   21.736455  18.287890
## 29   13      R70   12.7527374   14.743597  13.748167
## 30   13      R98   17.9034635   24.295281  21.099372
## 31   14      A10   16.0353867   20.806689  18.421038
## 32   14     B100   14.8758324   16.482712  15.679272
## 33   14     R102    9.5848380   17.897359  13.741099
## 34   14     R128   11.5441914   18.807803  15.175997
## 35   14     R133   24.4394863   26.317757  25.378622
## 36   14     R161   11.3351757   16.721924  14.028550
## 37   14     R187   15.9765377   22.156071  19.066305
## 38   14      R20    9.1279616   35.869895  22.498929
## 39   14      R70   12.8683843   21.043290  16.955837
## 40   14      R98   22.7445324   26.529747  24.637140
## 41   15      A10   20.7800983   23.377751  22.078924
## 42   15     B100   17.7777752   20.518164  19.147970
## 43   15     R102   11.0573174   17.922030  14.489674
## 44   15     R128   16.7501844   21.295836  19.023010
## 45   15     R133   27.0393189   33.948830  30.494074
## 46   15     R161   17.1755797   20.612348  18.893964
## 47   15     R187   15.8387658   22.954208  19.396487
## 48   15      R20   25.1788370   33.253535  29.216186
## 49   15      R70   18.9601854   22.863978  20.912081
## 50   15      R98   29.0498999   31.955600  30.502750
## 51   16      A10   20.0577561   25.806257  22.932007
## 52   16     B100   20.0413990   23.203604  21.622501
## 53   16     R102   14.8670912   21.144546  18.005818
## 54   16     R128   16.5582461   30.333804  23.446025
## 55   16     R133   33.0362095   35.768486  34.402348
## 56   16     R161   16.4543384   23.558591  20.006465
## 57   16     R187   22.3958994   28.061383  25.228641
## 58   16      R20    9.0056082   50.912573  29.959091
## 59   16      R70   20.2903256   31.959948  26.125137
## 60   16      R98   32.7722175   35.627625  34.199921
## 61   17      A10   27.4881214   32.177027  29.832574
## 62   17     B100   25.7375286   29.145299  27.441414
## 63   17     R102    9.3832005   41.668709  25.525955
## 64   17     R128   26.8195976   35.438632  31.129115
## 65   17     R133   35.0134141   35.046324  35.029869
## 66   17     R161   18.5118930   41.274711  29.893302
## 67   17     R187   26.7251684   31.553360  29.139264
## 68   17      R20   34.5958844   35.180457  34.888171
## 69   17      R70   27.8333370   32.092941  29.963139
## 70   17      R98   34.9761994   35.063296  35.019748
## 71   18      A10   25.4649891   34.113739  29.789364
## 72   18     B100   29.5995162   33.060621  31.330069
## 73   18     R102   23.9458634   35.780013  29.862938
## 74   18     R128   28.8943167   33.278893  31.086605
## 75   18     R161   22.7950766   33.274113  28.034595
## 76   18     R187   28.6642353   38.039275  33.351755
## 77   18      R20   34.8319033   35.217714  35.024809
## 78   18      R70   32.5635770   36.332210  34.447894
## 79   18      R98   35.0377193   35.072626  35.055172
## 80   19      A10   29.9565044   33.715470  31.835987
## 81   19     B100   28.7587854   32.088308  30.423547
## 82   19     R102   18.6004422   36.868409  27.734426
## 83   19     R128   35.0349633   35.062368  35.048666
## 84   19     R133   35.0369049   35.055222  35.046063
## 85   19     R161   27.5528428   35.523189  31.538016
## 86   19     R187   30.2853896   36.052235  33.168812
## 87   19      R20   35.0431394   35.064445  35.053792
## 88   19      R70   32.8273407   35.931809  34.379575
## 89   20      A10   30.7762330   35.328096  33.052164
## 90   20     B100   31.4095265   35.224147  33.316837
## 91   20     R102   29.8675676   36.244858  33.056213
## 92   20     R128   31.8182665   36.308301  34.063284
## 93   20     R133   35.0292025   35.062924  35.046063
## 94   20     R161   27.6114760   34.513152  31.062314
## 95   20     R187   33.7013663   35.619516  34.660441
## 96   20      R20   35.0065633   35.093660  35.050112
## 97   20      R70   31.1659064   36.873618  34.019762
## 98   20      R98   35.0336568   35.066567  35.050112
## 99   21      A10   34.7350643   43.619676  39.177370
## 100  21     B100   33.9185376   40.511408  37.214973
## 101  21     R102   20.4335385   53.700974  37.067256
## 102  21     R128   45.7878729   59.113583  52.450728
## 103  21     R133   69.5748771   75.552885  72.563881
## 104  21     R161   30.8383653   46.826065  38.832215
## 105  21     R187   41.2621776   55.259613  48.260895
## 106  21      R20   56.3428640   68.941224  62.642044
## 107  21      R70   45.2148076   52.242602  48.728705
## 108  21      R98   56.5053136   69.874481  63.189898
## 109  22      A10   42.1333456   51.279316  46.706331
## 110  22     B100   35.1894676   46.590960  40.890214
## 111  22     R102   32.4270279   54.660958  43.543993
## 112  22     R128   42.7607416   57.866538  50.313640
## 113  22     R133   64.7341342   95.014087  79.874111
## 114  22     R161   31.9080496   39.897530  35.902790
## 115  22     R187   49.2083816   68.125731  58.667056
## 116  22      R20   45.5047607  113.876490  79.690625
## 117  22      R70   37.0178551   70.222627  53.620241
## 118  22      R98   64.6229237   77.916294  71.269609
## 119  23      A10   46.1205983   57.372009  51.746303
## 120  23     B100   39.0566070   49.213411  44.135009
## 121  23     R102   22.1355965   62.113193  42.124395
## 122  23     R128   56.6189532   68.402234  62.510593
## 123  23     R133   80.1898844   87.496497  83.843190
## 124  23     R161   39.5976979   48.909886  44.253792
## 125  23     R187   44.9050613   74.178923  59.541992
## 126  23      R20   65.9802051   85.598085  75.789145
## 127  23      R70   58.6470448   67.651642  63.149343
## 128  23      R98   60.6980969   79.556520  70.127309
## 129  25      A10   55.9173941   66.531485  61.224439
## 130  25     B100   42.1436759   60.516430  51.330053
## 131  25     R102   37.5833391   74.825638  56.204488
## 132  25     R128   59.0014536   77.608288  68.304871
## 133  25     R133   91.5153533  106.509863  99.012608
## 134  25     R161   41.7599259   51.527936  46.643931
## 135  25     R187   69.4178315   80.320776  74.869304
## 136  25      R20   47.3113464  137.062553  92.186950
## 137  25      R70   54.8053892   93.180449  73.992919
## 138  25      R98   79.0174138   95.332386  87.174900
## 139  26      A10   60.7309054   74.484346  67.607626
## 140  26     B100   51.9573893   67.031311  59.494350
## 141  26     R102   29.0283221   90.704621  59.866472
## 142  26     R128   69.7432021   85.949025  77.846114
## 143  26     R133  102.1131687  102.209599 102.161384
## 144  26     R161   44.7659963   77.522223  61.144110
## 145  26     R187   64.2493428   83.643787  73.946565
## 146  26      R20   95.4818290  103.551591  99.516710
## 147  26      R70   80.8950037   90.615600  85.755302
## 148  26      R98   67.3596660   96.850149  82.104907
## 149  27      A10   68.1673027   79.966027  74.066665
## 150  27     B100   52.3930899   73.531038  62.962064
## 151  27     R102   46.4307081   91.299260  68.864984
## 152  27     R128   69.7163940   92.330224  81.023309
## 153  27     R133  102.0534520  102.623594 102.338523
## 154  27     R161   47.0038080   59.572263  53.288036
## 155  27     R187   79.9731196   87.133834  83.553477
## 156  27      R20   94.6611877  107.250057 100.955622
## 157  27      R70   67.7411474  102.378830  85.059989
## 158  27      R98   92.2825030  102.888066  97.585284
## 159  28      A10   72.4019654   84.816088  78.609027
## 160  28     B100   60.2668834   76.201220  68.234052
## 161  28     R102   38.3458163  100.670481  69.508149
## 162  28     R128   88.7635776   97.020262  92.891920
## 163  28     R133  102.0427075  102.412501 102.227604
## 164  28     R161   54.3065643   86.175685  70.241125
## 165  28     R187   68.2614296   91.127546  79.694488
## 166  28      R20  101.9335479  102.697144 102.315346
## 167  28      R70   96.9145446  101.216624  99.065584
## 168  28      R98   81.5682665  106.133203  93.850735
## 169  29      A10   77.9931551   87.733707  82.863431
## 170  29     B100   64.8732783   90.167093  77.520186
## 171  29     R102   51.4886408  103.817301  77.652971
## 172  29     R128   78.2635038  103.522285  90.892894
## 173  29     R133  102.1483036  102.853222 102.500763
## 174  29     R161   60.8247415   72.598147  66.711444
## 175  29     R187   85.3621398   92.131055  88.746598
## 176  29      R20  102.1549445  102.730696 102.442820
## 177  29      R70   81.5101446  109.747356  95.628750
## 178  29      R98  100.9210632  102.477379 101.699221
## 179  30      A10   79.7975021   93.165534  86.481518
## 180  30     B100   69.9817638   86.831826  78.406795
## 181  30     R102   52.6392582  105.567685  79.103472
## 182  30     R128  100.2996586  102.191014 101.245336
## 183  30     R133  102.1362130  102.517656 102.326935
## 184  30     R161   72.2477821   89.354608  80.801195
## 185  30     R187   75.6166966   95.754845  85.685771
## 186  30      R20  102.3143461  102.675591 102.494969
## 187  30      R70  102.0282615  102.642163 102.335212
## 188  30      R98   93.7598548  105.298398  99.529126
## 189  31      A10   87.1296524   94.989746  91.059699
## 190  31     B100   76.9270063   98.740991  87.833999
## 191  31     R102   62.8067625  104.995504  83.901133
## 192  31     R128   89.1628356  105.911162  97.536999
## 193  31     R133  102.3840680  102.756520 102.570294
## 194  31     R161   68.6018081   80.383543  74.492676
## 195  31     R187   89.8228715   96.405196  93.114034
## 196  31      R20  101.8550011  103.378296 102.616648
## 197  31      R70   87.6302155  109.560625  98.595420
## 198  31      R98  102.2831228  102.641146 102.462134
## 199  32      A10   86.6320427  101.920801  94.276422
## 200  32     B100   82.8496624   97.181019  90.015341
## 201  32     R102   66.9055230  110.225409  88.565466
## 202  32     R128  105.0325661  107.488859 106.260713
## 203  32     R133  107.1520311  107.382340 107.267186
## 204  32     R161   77.4568438   96.766731  87.111787
## 205  32     R187   85.8563435   97.724316  91.790330
## 206  32      R20  107.2913801  107.379954 107.335667
## 207  32      R70  107.0598209  107.278888 107.169355
## 208  32      R98  105.1210435  108.008120 106.564582
## 209  33      A10   96.0947018  103.735921  99.915311
## 210  33     B100   89.2731829  105.354814  97.313999
## 211  33     R102   75.8742679  110.271342  93.072805
## 212  33     R128   96.1577692  110.345571 103.251670
## 213  33     R133  107.1725115  107.361860 107.267186
## 214  33     R161   81.8858866   88.986237  85.436062
## 215  33     R187   93.8716625   99.386701  96.629182
## 216  33      R20  107.0615755  107.597308 107.329442
## 217  33      R70   92.3612921  114.279011 103.320152
## 218  33      R98  107.1300106  107.383609 107.256810

Treatment differences in height

Statistical analysis for height differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.

analyze_height = function(vis.data, genotype) {
  days = c()
  diff.low = c()
  diff.up = c()
  pvals = c()
  for(day in levels(factor(as.integer(vis.data$dap)))) {
    day = as.integer(day)
    control = vis.data[(as.integer(vis.data$dap) == day |
                        as.integer(vis.data$dap) == day + 1) &
                        vis.data$genotype == genotype &
                        vis.data$treatment == 100,]$height_above_bound
    drought = vis.data[(as.integer(vis.data$dap) == day |
                        as.integer(vis.data$dap) == day + 1) &
                        vis.data$genotype == genotype &
                        vis.data$treatment == 33,]$height_above_bound
    test = t.test(x=control, y=drought)
    days = c(days, day)
    diff.low = c(diff.low, test$conf.int[1])
    diff.up = c(diff.up, test$conf.int[2])
    pvals = c(pvals, test$p.value)
  }
  results = data.frame(dap=as.numeric(days),
                       conf.int.low=diff.low,
                       conf.int.up=diff.up,
                       pvalue=pvals)
  return(results)
}

Test for treatment effect on two-day intervals

a10.height.results = analyze_height(vis.data, 'A10')
b100.height.results = analyze_height(vis.data, 'B100')

Control for multiple testing by controlling the FDR

qvalues.height = qvalue(c(a10.height.results$pvalue, b100.height.results$pvalue))
a10.height.results$qvalue = qvalues.height$qvalues[1:nrow(a10.height.results)]
b100.height.results$qvalue = qvalues.height$qvalues[
  (nrow(a10.height.results)+1):(nrow(a10.height.results) + nrow(b100.height.results))]

Assign genotypes for merged table

a10.height.results$genotype = 'A10'
b100.height.results$genotype = 'B100'
print(a10.height.results)
##    dap conf.int.low conf.int.up       pvalue       qvalue genotype
## 1   11    1.3757087    5.139721 1.622327e-03 1.057500e-03      A10
## 2   12   -0.7380484    4.515406 1.490756e-01 5.492423e-02      A10
## 3   13   -0.5912495    3.543282 1.562446e-01 5.516696e-02      A10
## 4   14   -1.1319577    4.329657 2.381965e-01 6.960212e-02      A10
## 5   15   -0.1299653    5.167933 6.143652e-02 3.062407e-02      A10
## 6   16   -1.3607720    7.080051 1.715524e-01 5.519138e-02      A10
## 7   17   -2.9248819    5.400151 5.353411e-01 1.191457e-01      A10
## 8   18   -0.4959663    6.673738 8.769306e-02 3.911081e-02      A10
## 9   19    1.4409925    7.969806 6.845578e-03 4.143499e-03      A10
## 10  20    3.6149867   12.029236 6.790664e-04 4.795304e-04      A10
## 11  21    6.7803624   17.104831 5.343934e-05 6.694424e-05      A10
## 12  22   10.2910019   21.102986 1.471517e-06 1.246954e-05      A10
## 13  23    9.2657053   24.761796 2.155551e-04 1.660545e-04      A10
## 14  25    7.2790852   20.484417 1.699507e-04 1.600168e-04      A10
## 15  26    8.5545706   23.844353 1.911306e-04 1.619628e-04      A10
## 16  27   11.5264141   26.176502 1.946187e-05 3.758465e-05      A10
## 17  28   10.7457158   25.976630 5.606520e-05 6.694424e-05      A10
## 18  29   11.4473342   24.918441 7.632409e-06 3.233827e-05      A10
## 19  30   11.9787099   27.114832 2.217662e-05 3.758465e-05      A10
## 20  31   12.9567838   28.606319 1.247030e-05 3.522418e-05      A10
## 21  32   12.2483513   29.889421 6.320013e-05 6.694424e-05      A10
## 22  33    3.4267550   38.584764 2.690657e-02 1.518525e-02      A10
print(b100.height.results)
##    dap conf.int.low conf.int.up     pvalue     qvalue genotype
## 1   11   -2.1888257   0.4486175 0.18236609 0.05519138     B100
## 2   12   -2.2019131   1.9489490 0.89939162 0.17724157     B100
## 3   13   -3.7808622   0.5052064 0.12627973 0.05095650     B100
## 4   14   -2.7496601   2.5699624 0.94357005 0.18172165     B100
## 5   15   -4.8854099   0.8999123 0.16277335 0.05517324     B100
## 6   16   -2.8974999   3.9598806 0.74490405 0.15029213     B100
## 7   17   -5.0026064   2.9536026 0.58132070 0.12014816     B100
## 8   18   -1.9644558   5.1388196 0.35534030 0.09713326     B100
## 9   19   -2.3699180   4.9749886 0.46324842 0.10904271     B100
## 10  20   -0.8187179   6.7951985 0.11836288 0.05014997     B100
## 11  21   -8.9817009   4.6032442 0.50430459 0.11549851     B100
## 12  22   -8.3643015   4.7370840 0.57202885 0.12014816     B100
## 13  23  -11.0131272   5.4448269 0.46303807 0.10904271     B100
## 14  25  -13.4045213   2.5892801 0.17774764 0.05519138     B100
## 15  26  -17.3140056  -1.0219101 0.02867192 0.01518525     B100
## 16  27  -14.3161656   2.1808693 0.14381616 0.05492423     B100
## 17  28  -17.0011012   1.2317531 0.08754828 0.03911081     B100
## 18  29  -11.9653330   5.2377910 0.43167803 0.10904271     B100
## 19  30  -12.3111302   4.5263257 0.35381611 0.09713326     B100
## 20  31   -4.5854432  10.1681899 0.44686134 0.10904271     B100
## 21  32   -5.5221448  10.1865601 0.54834969 0.11914565     B100
## 22  33   -6.1404245  14.8562724 0.38830762 0.10282794     B100

Output the A10 and B100 tables to the same file

write.table(a10.height.results, file='height.stats.csv', sep = ',',
            row.names = FALSE, append = FALSE)
write.table(b100.height.results, file='height.stats.csv', sep = ',',
            row.names = FALSE, append = TRUE, col.names = FALSE)

Biomass analysis

Plot fresh-weight biomass for S. viridis and S. italica water treatments 100% and 33% full-capacity.

biomass.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
                                vis.data$genotype == 'B100') &
                               (vis.data$treatment == 100 |
                                vis.data$treatment == 33),],
                     aes(x=dap, y=fw_biomass, color=factor(group))) +
                     geom_point(size=2.5) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(-1,41),
                                        name="Estimated biomass (g)") +
                     theme_bw() +
                     theme(legend.position=c(0.25,0.75),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Genotype-treatment")
biomass.plot = biomass.plot + labs(title='Figure 4B')
print(biomass.plot)

Plot fresh-weight biomass for all genotypes.

biomass.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
                                    vis.data$treatment == 33,],
                           aes(x=dap, y=fw_biomass, color=factor(treatment))) +
                     geom_point(size=1) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(-1,41),
                                        name="Estimated biomass (g)") +
                     theme_bw() +
                     theme(legend.position=c(0.8,0.1),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Treatment") +
                     facet_wrap(~ genotype, ncol = 3)
biomass.facet.plot = biomass.facet.plot + labs(title='Supplemental Figure S3')
print(biomass.facet.plot)

Treatment differences in biomass

Statistical analysis for fresh-weight biomass differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.

analyze_biomass = function(vis.data, genotype) {
  days = c()
  diff.low = c()
  diff.up = c()
  pvals = c()
  for(day in levels(factor(as.integer(vis.data$dap)))) {
    day = as.integer(day)
    control = vis.data[(as.integer(vis.data$dap) == day |
                        as.integer(vis.data$dap) == day + 1) &
                        vis.data$genotype == genotype &
                        vis.data$treatment == 100,]$fw_biomass
    drought = vis.data[(as.integer(vis.data$dap) == day |
                        as.integer(vis.data$dap) == day + 1) &
                        vis.data$genotype == genotype &
                        vis.data$treatment == 33,]$fw_biomass
    # If the control group has a lot more replicates (e.g. A10 and B100), randomly sample to drought sample size
    if (genotype == 'A10' | genotype == 'B100') {
      control = sample(control, size = length(drought))
    }
    test = t.test(x=control, y=drought)
    days = c(days, day)
    diff.low = c(diff.low, test$conf.int[1])
    diff.up = c(diff.up, test$conf.int[2])
    pvals = c(pvals, test$p.value)
  }
  results = data.frame(dap=as.numeric(days),
                       conf.int.low=diff.low,
                       conf.int.up=diff.up,
                       pvalue=pvals)
  return(results)
}

Test for treatment effect on two-day intervals

a10.biomass.results = analyze_biomass(vis.data, 'A10')
b100.biomass.results = analyze_biomass(vis.data, 'B100')

Control for multiple testing by controlling the FDR

qvalues.biomass = qvalue(c(a10.biomass.results$pvalue, b100.biomass.results$pvalue))
a10.biomass.results$qvalue = qvalues.biomass$qvalues[1:nrow(a10.biomass.results)]
b100.biomass.results$qvalue = qvalues.biomass$qvalues[
  (nrow(a10.biomass.results)+1):(nrow(a10.biomass.results) + 
                                 nrow(b100.biomass.results))]

Assign genotypes for merged table

a10.biomass.results$genotype = 'A10'
b100.biomass.results$genotype = 'B100'
print(a10.biomass.results)
##    dap conf.int.low conf.int.up       pvalue       qvalue genotype
## 1   11   0.02914367   0.2201204 1.323444e-02 1.515363e-02      A10
## 2   12  -0.14788575   0.1756949 8.610677e-01 5.049914e-01      A10
## 3   13  -0.16165563   0.2625467 6.306207e-01 3.790869e-01      A10
## 4   14  -0.18405590   0.4646706 3.796492e-01 2.552294e-01      A10
## 5   15  -0.27354779   0.6885945 3.821226e-01 2.552294e-01      A10
## 6   16  -0.39537050   0.6661931 5.981374e-01 3.784843e-01      A10
## 7   17   0.10518290   3.1083996 3.780833e-02 3.952667e-02      A10
## 8   18  -0.13454475   2.2779082 7.835338e-02 7.138478e-02      A10
## 9   19   0.55639946   3.2418022 7.520995e-03 9.518142e-03      A10
## 10  20   1.79793954   4.9902898 2.723816e-04 4.093438e-04      A10
## 11  21   2.01188934   6.2302349 6.230113e-04 8.812062e-04      A10
## 12  22   4.48502027   8.5942426 2.003136e-06 4.013837e-06      A10
## 13  23   4.38451816  10.3131422 2.161645e-04 3.465163e-04      A10
## 14  25   6.34140762  10.9450081 1.001322e-07 3.009637e-07      A10
## 15  26   6.34064811  11.1993657 2.033517e-07 5.432952e-07      A10
## 16  27   8.43389988  12.1214602 4.849913e-11 2.332354e-10      A10
## 17  28   9.48931032  12.8236353 6.467846e-13 7.776071e-12      A10
## 18  29   8.43000027  11.5285952 1.588593e-12 1.041987e-11      A10
## 19  30  10.37389366  13.1633241 1.396552e-14 3.358054e-13      A10
## 20  31   9.37549132  12.8536239 1.733372e-12 1.041987e-11      A10
## 21  32   8.77790792  12.5759125 8.718795e-11 3.494103e-10      A10
## 22  33   9.13675701  14.3597517 1.794587e-06 4.013837e-06      A10
print(b100.biomass.results)
##    dap conf.int.low conf.int.up       pvalue       qvalue genotype
## 1   11   -0.1298704 -0.02014007 9.958668e-03 1.197297e-02     B100
## 2   12   -0.0805700  0.14440699 5.635119e-01 3.662114e-01     B100
## 3   13   -0.2655002  0.05797019 1.994123e-01 1.453010e-01     B100
## 4   14   -0.2982824  0.27327839 9.285308e-01 5.315909e-01     B100
## 5   15   -0.5579089  0.03474439 8.015652e-02 7.138478e-02     B100
## 6   16   -0.4525040  0.74118118 6.175631e-01 3.790869e-01     B100
## 7   17   -1.0828481  1.08346063 9.995283e-01 5.462267e-01     B100
## 8   18   -0.2717842  1.41231582 1.737604e-01 1.347782e-01     B100
## 9   19   -1.1177500  1.15611705 9.725059e-01 5.438189e-01     B100
## 10  20   -0.2312011  2.65901500 9.508123e-02 8.165209e-02     B100
## 11  21   -0.2643127  2.63411000 1.034254e-01 8.575509e-02     B100
## 12  22   -0.5396748  3.20099463 1.535486e-01 1.230708e-01     B100
## 13  23   -2.9434030  7.26769789 3.233459e-01 2.286751e-01     B100
## 14  25   -1.2656478  5.61185394 1.982686e-01 1.453010e-01     B100
## 15  26   -0.1065401  5.98201981 5.765873e-02 5.545690e-02     B100
## 16  27    0.1611022  6.85503580 4.110506e-02 4.118268e-02     B100
## 17  28    0.7754796  9.21814913 2.356266e-02 2.575326e-02     B100
## 18  29    4.0254297 10.07486231 1.391536e-04 2.389995e-04     B100
## 19  30    6.9552625 11.09939185 7.806236e-09 2.681477e-08     B100
## 20  31    5.5617792 10.73510877 2.262272e-06 4.184388e-06     B100
## 21  32    3.0137106 11.18548970 2.080233e-03 2.778882e-03     B100
## 22  33    8.1298505 13.84599687 1.865143e-06 4.013837e-06     B100

Output the A10 and B100 tables to the same file

write.table(a10.biomass.results, file='biomass.stats.csv', sep = ',',
            row.names = FALSE, append = FALSE)
write.table(b100.biomass.results, file='biomass.stats.csv', sep = ',',
            row.names = FALSE, append = TRUE, col.names = FALSE)

Biomass response to drought

Calculate the difference in biomass per genotype per day between the 33% and 100% water treatment groups.

drought_response = function(vis.data, genotypes) {
  # Initialize drought response data frame with days
  drought_resp = data.frame(day=c(
    min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))))
  
  # Initialize genotypes
  for(g in genotypes) {
    control = paste(g,'control',sep='')
    drought = paste(g,'drought',sep='')
    
    # Calculate median biomass per treatment per day
    drought_resp[,paste(control)] = 0
    drought_resp[,paste(drought)] = 0
    for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
      drought_resp[drought_resp$day == d, paste(control)] = 
        median(vis.data[vis.data$genotype == g &
                        vis.data$treatment == 100 &
                        as.integer(vis.data$dap) == d,]$fw_biomass)
      drought_resp[drought_resp$day == d, paste(drought)] = 
        median(vis.data[vis.data$genotype == g &
                        vis.data$treatment == 33 &
                        as.integer(vis.data$dap) == d,]$fw_biomass)
    }
  }
  drought_resp = drought_resp[-14,]
  
  # Calculate biomass loss to drought
  days = c()
  response = c()
  genotype=c()
  for(r in 1:nrow(drought_resp)) {
    days = c(days, rep(drought_resp[r,]$day,length(genotypes)))
    for(g in genotypes) {
      control = paste(g,'control',sep='')
      drought = paste(g,'drought',sep='')
      response = c(response, drought_resp[r,paste(drought)] -
                   drought_resp[r,paste(control)])
      genotype = c(genotype, g)
    }
  }
  dr.df = data.frame(day=days,response=response,genotype=as.factor(genotype))
  return(dr.df)
}

A10 vs B100 drought responses

drought.response.parents = drought_response(vis.data,c('A10','B100'))

Plot drought response results

genotype.colors = c('#E6A024','#4AB859')
resp.plot = ggplot(drought.response.parents,
                   aes(x=day, y=response, group=genotype, color=genotype)) +
                   geom_point(size=2.5) +
                   geom_smooth(method="loess") +
                   scale_x_continuous(name="Days after planting") +
                   scale_y_continuous(name="Reduced accumulated biomass (g)") +
                   theme_bw() +
                   theme(legend.position=c(0.8,0.8),
                         axis.title.x=element_text(face="bold"),
                         axis.title.y=element_text(face="bold")) +
                   scale_colour_manual(values=genotype.colors) +
                   labs(color='Genotype')
resp.plot = resp.plot + labs(title='Figure 4C')
print(resp.plot)

Growth curve analysis

In this experiment plants appear to follow relatively basic asymptotic growth patterns. Here we use non-linear least squares regression analysis to model a three-component logistic growth function for each genotype x treatment group.

First load three functions for non-linear growth curve analysis reproduced from:

Paine CET, Marthews TR, Vogt DR, Purves D, Rees M, Hector A, Turnbull LA (2012) How to fit nonlinear plant growth models and calculate growth rates: an update for ecologists. Methods in Ecology and Evolution 3: 245–256. doi: 10.1111/j.2041-210X.2011.00155.x.

output.logis.nlsList <- function(fit, times, CI = F, LOG = F, alpha = 0.05){
  coef <- coef(fit)
  params <- transform_param.logis(coef)
  rates <- list()
  groups <- rownames(params)
  n.groups <- nrow(coef)
  
  # compute rates for each group seperately
  rates <- list()
  for(i in 1:(n.groups)){
    K <- params[i,1]; r <- params[i,2]; M0 <- params[i,3]
    rates[[i]] = data.frame(
      times = times,
      M    = (M0*K)/(M0+(K-M0)*exp(-r*times)),
      AGR  = (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
    )
    rates[[i]]$RGRt <- rates[[i]]$AGR/rates[[i]]$M
    rates[[i]]$RGRm <- r*(1 - rates[[i]]$M/K)
    if(LOG == T){
      rates[[i]]$RGRt <- rates[[i]]$AGR
      rates[[i]]$RGRm <- r*rates[[i]]$M*(1-rates[[i]]$M/K)
      rates[[i]]$AGR  <- rates[[i]]$AGR*exp(rates[[i]]$M)
    }
    # commute CIs for each group's estaimates, if desired
    if(CI == T){
      cov   <- summary(fit)$cov[i,,]
      x <- y <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=cov))
      x$K  <- y[,1]
      x$r  <- 1/y[,3]
      x$M0 <- y[,1]/(1 + exp(y[,2]/y[,3])) 
      M <- AGR <- RGRt <- RGRm <- matrix(NA, ncol = length(times), nrow = nrow(x))
      for(j in 1:nrow(x)){
        K <- x[j,4]; r <- x[j,5]; M0 <- x[j,6]
        M[j,]     <- (M0*K)/(M0+(K-M0)*exp(-r*times))
        AGR[j,]   <- (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
        RGRt[j,]  <- AGR[j,]/M[j,]
        RGRm[j,]  <- r*(1 - M[j,]/K)
        if(LOG ==T){
          RGRt[j,] <- AGR[j,]
          RGRm[j,] <- r*M[j,]*(1 - M[j,]/K)
          AGR[j,]  <- AGR[j,]*exp(M[j,])
        }
      }
    }
    CIs <- summarizer(list(M = M, AGR = AGR, RGRt = RGRt, RGRm = RGRm), alpha)
    rates[[i]]  <-  cbind(rates[[i]], CIs)
  }
  names(rates) <- rownames(params)
  
  # now compute differences among groups
  diffs <- list()
  for(i in 1:(n.groups-1)){
    Ki <- params[i,1]; ri <- params[i,2]; M0i <- params[i,3]
    for(j in (i+1):n.groups){
      Kj <- params[j,1]; rj <- params[j,2]; M0j <- params[j,3]
      diffs.ij = data.frame(
        times = times,
        diffM     = rates[[i]]$M    - rates[[j]]$M,
        diffAGR   = rates[[i]]$AGR  - rates[[j]]$AGR,
        diffRGRt  = rates[[i]]$RGRt - rates[[j]]$RGRt
      )
      # comparing RGRm has to be done on a biomass basis. So it needs special treatment. First, we aev to know what range of biomasses are shared between groups.
      Mmin <- max(min(rates[[i]]$M), min(rates[[j]]$M)) # yieds the range of overlapping masses between groups i & j
      Mmax <- min(max(rates[[i]]$M), max(rates[[j]]$M))
      diffs.ij$Mseq <- seq(Mmin, Mmax, length = 100)
      if(LOG == F){
        diffs.ij$diffRGRm  <- ri*(1 - diffs.ij$Mseq/Ki) - rj*(1 - diffs.ij$Mseq/Kj)
      } else{
        diffs.ij$diffRGRm <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki) - rj*Mseq*(1 - diffs.ij$Mseq/Kj)
      }    
      
    }
    if(CI == T){
      # get params for group i
      covi   <- summary(fit)$cov[i,,]
      xi <- yi <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=covi))
      xi$K  <- yi[,1]
      xi$r  <- 1/yi[,3]
      xi$M0 <- yi[,1]/(1 + exp(yi[,2]/yi[,3])) 
      
      # get params for group j
      covj   <- summary(fit)$cov[j,,]
      xj <- yj <- data.frame(rmvnorm(n=1000, mean=c(coef[j, 1], coef[j, 2], coef[j, 3]), sigma=covj))
      xj$K  <- yj[,1]
      xj$r  <- 1/yj[,3]
      xj$M0 <- yj[,1]/(1 + exp(yj[,2]/yj[,3])) 
      
      # now compute diffs for each random set of drawn parameters
      Mi <- Mj <- AGRi <- AGRj <- RGRti <- RGRtj <- RGRmi <- RGRmj <- diffM <- diffAGR <- diffRGRt <- diffRGRm <- matrix(NA, ncol = length(times), nrow = nrow(xi))
      for(k in 1:nrow(xi)){
        Ki <- xi[k,4]; ri <- xi[k,5]; M0i <- xi[k,6]
        Kj <- xj[k,4]; rj <- xj[k,5]; M0j <- xj[k,6]
        Mi[k,]    <- (M0i*Ki)/(M0i+(Ki-M0i)*exp(-ri*times))
        Mj[k,]    <- (M0j*Kj)/(M0j+(Kj-M0j)*exp(-rj*times))
        AGRi[k,]  <- (ri*M0i*Ki*(Ki-M0i)*exp(-ri*times))/(M0i+(Ki-M0i)*exp(-ri*times))^2
        AGRj[k,]  <- (rj*M0j*Kj*(Kj-M0j)*exp(-rj*times))/(M0j+(Kj-M0j)*exp(-rj*times))^2
        RGRti[k,] <- AGRi[k,]/Mi[k,]
        RGRtj[k,] <- AGRj[k,]/Mj[k,]
        RGRmi[k,] <- ri*(1 - diffs.ij$Mseq/Ki)
        RGRmj[k,] <- rj*(1 - diffs.ij$Mseq/Kj)
        if(LOG == T){
          RGRti[k,] <- AGRi[k,]
          RGRtj[k,] <- AGRj[k,]
          RGRmi[k,] <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki)
          RGRmj[k,] <- rj*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Kj)
          AGRi[k,]  <- AGRi[k,]*exp(Mi[k,])
          AGRj[k,]  <- AGRj[k,]*exp(Mj[k,])
        }
        diffM[k,]    <- Mi[k,]    - Mj[k,]
        diffAGR[k,]  <- AGRi[k,]  - AGRj[k,]
        diffRGRt[k,] <- RGRti[k,] - RGRtj[k,]
        diffRGRm[k,] <- RGRmi[k,] - RGRmj[k,]
      }
      CIs <- summarizer(list(diffM = diffM, diffAGR = diffAGR, diffRGRt = diffRGRt, diffRGRm = diffRGRm), alpha)
      diffs[[paste(groups[i], groups[j], sep = "_")]]  <-  cbind(diffs.ij, CIs)
    } else{
      diffs[[paste(groups[i], groups[j], sep = "_")]]  <-  diffs.ij
    }
  } # end loop over pairwise combinations of groups
  
  out <- list(params = params, rates = rates, diffs = diffs)
  return(out)
}

transform_param.logis <- function(coef){
  K = coef[1]
  r = 1/(coef[3])
  M0 =  K/(1 + exp(coef[2]/coef[3])) #untransform best-fit parameters to K, r and M0
  if(is.data.frame(K)){
    out <- cbind(K, r, M0)
  } else {
    out <- c(K, r, M0)
  }
  names(out) <- c("K", "r", "M0")
  return(out)
}

# this function returns confidence envelopes around growth trajectories, and growth rates. 
summarizer <- function(dat, alpha){
  n <- length(dat)
  quantiles <- c(alpha/2, 1-(alpha/2))
  CIs <- data.frame(matrix(NA, ncol(dat[[1]]), n*2))
  names(CIs) <- paste(rep(names(dat), each = 2), c("lo", "hi"), sep = ".")
  for(i in 1:n){
    CIs[,(2*i-1):(2*i)] <- t(apply(dat[[i]],    2, quantile, quantiles, na.rm = T))
  }
  return(CIs)
}

Subset the VIS data

sub.vis.data = vis.data[(vis.data$treatment == 100 | vis.data$treatment == 33),
                        c("plant_id","dap","fw_biomass","group")]
sub.vis.data$group = factor(sub.vis.data$group)

Group data

grouped.sub.vis.data = groupedData(fw_biomass ~ dap | group, sub.vis.data)

Fit three-component logistic functions for each genotype x treatment group

fit.logis = nlsList(fw_biomass ~ SSlogis(dap, Asym, xmid, scal),
                    data = grouped.sub.vis.data)

Output estimated growth rate parameters at even time intervals

est.interval = seq(min(sub.vis.data$dap), max(sub.vis.data$dap), length = 100)
out.fit.logis = output.logis.nlsList(fit.logis, times = est.interval,
                                     CI = TRUE, LOG = FALSE, alpha = 0.05)

Compute the time and magnitude of maximum growth rate

logis.results = data.frame(group=levels(sub.vis.data$group))
logis.results$max.AGR = 0
logis.results$max.AGR.dap = 0
group.order = names(out.fit.logis$rates)
for(i in 1:length(group.order)) {
  logis.results[logis.results == group.order[i],]$max.AGR.dap = 
    out.fit.logis$rates[[i]]$times[out.fit.logis$rates[[i]]$AGR == 
                                             max(out.fit.logis$rates[[i]]$AGR)]
  logis.results[logis.results == group.order[i],]$max.AGR =
    max(out.fit.logis$rates[[i]]$AGR)
}
print(logis.results)
##       group  max.AGR max.AGR.dap
## 1   A10-100 2.994229    24.23556
## 2    A10-33 1.737080    25.36054
## 3  B100-100 2.389559    26.26053
## 4   B100-33 1.631987    24.68555
## 5  R102-100 2.593812    28.06050
## 6   R102-33 1.826986    26.48553
## 7  R128-100 2.995409    25.36054
## 8   R128-33 1.781943    24.23556
## 9  R133-100 3.267807    23.11057
## 10  R133-33 2.207510    23.33557
## 11 R161-100 2.858742    25.58554
## 12  R161-33 1.750648    24.68555
## 13 R187-100 2.444363    24.46055
## 14  R187-33 1.617136    24.46055
## 15  R20-100 3.240594    22.88557
## 16   R20-33 1.869087    22.66058
## 17  R70-100 2.878034    23.78556
## 18   R70-33 1.592505    23.33557
## 19  R98-100 2.509740    24.23556
## 20   R98-33 1.481989    24.68555

Plot logistic growth models for S. viridis and S. italica and absolute growth rates over time

parent.groups = c("A10-100", "A10-33", "B100-100", "B100-33")
group.colors = c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")

# Biomass over time
gca.plot = ggplot(vis.data[(vis.data$genotype == 'A10' | vis.data$genotype == 'B100') &
                 (vis.data$treatment == 100 | vis.data$treatment == 33),],
                  aes(x=dap, y=fw_biomass, color=factor(group))) +
                  geom_point(size=2.5) +
                  scale_x_continuous(name="Days after planting") +
                  scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
                  theme_bw() +
                  theme(legend.position=c(0.2,0.8),
                        axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold")) +
                  scale_colour_manual(name='Genotype-treatment',
                                      values=group.colors, labels=parent.groups)

agr.plot = ggplot() + 
           scale_x_continuous(name="Days after planting") +
           scale_y_continuous(name="Absolute growth rate (g/day)") +
           theme_bw() +
           theme(legend.position=c(0.2,0.8),
                        axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold")) +
           scale_colour_manual(name='Genotype-treatment',
                                      values=group.colors, labels=parent.groups)

# Add logistic growth models and confidence intervals
for(g in 1:length(parent.groups)) {
  group = parent.groups[g]
  i = match(group, group.order)
  group.rates = as.data.frame(out.fit.logis$rates[i])
  col.name = gsub('-', '.', group)
  conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
                            rev(group.rates[,paste(col.name, '.times', sep='')])),
                        y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
                            rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
  agr.conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
                            rev(group.rates[,paste(col.name, '.times', sep='')])),
                        y=c(group.rates[,paste(col.name, '.AGR.lo', sep='')],
                            rev(group.rates[,paste(col.name, '.AGR.hi', sep='')])))
  gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
                                     fill='gray60', color=NA, alpha=0.4) +
    geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
                                           y=paste(col.name, '.M', sep='')),
              color=group.colors[g])
  
  agr.plot = agr.plot +
    geom_polygon(data=agr.conf.int, aes(x=x, y=y),
                 fill='gray60', color=NA, alpha=0.4) +
    geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
                                           y=paste(col.name, '.AGR', sep='')),
              color=group.colors[g])
}
gca.plot = gca.plot + labs(title='Figure 4B')
print(gca.plot)

agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)

Plot logistic growth models for all Setaria genotypes

plot_gca = function(genotype) {
  groups = c(paste(genotype,'-100',sep=''), paste(genotype,'-33',sep=''))
  group.colors = c("#00BFC4", "#F8766D")
  # Biomass over time
  gca.plot = ggplot(vis.data[(vis.data$genotype == genotype) &
                   (vis.data$treatment == 100 | vis.data$treatment == 33),],
                    aes(x=dap, y=fw_biomass, color=factor(treatment))) +
                    geom_point(size=1) +
                    scale_x_continuous(name="Days after planting") +
                    scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
                    theme_bw() +
                    theme(legend.position='none',
                          axis.title.x=element_text(face="bold"),
                          axis.title.y=element_text(face="bold")) +
                    labs(color="Treatment") +
                    facet_grid(. ~ genotype)
  # Add logistic growth models and confidence intervals
  for(g in 1:length(groups)) {
    group = groups[g]
    i = match(group, group.order)
    group.rates = as.data.frame(out.fit.logis$rates[i])
    col.name = gsub('-', '.', group)
    conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
                            rev(group.rates[,paste(col.name, '.times', sep='')])),
                        y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
                            rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
    gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
                                     fill='gray60', color=NA, alpha=0.4) +
    geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
                                           y=paste(col.name, '.M', sep='')),
              color=group.colors[g])
  }
  return(gca.plot)
}

Multiple plot function modified from the Cookbook for R by Winston Chang: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/

multiplot <- function(..., plotlist=NULL, cols=1, layout=NULL) {
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols),byrow=TRUE)
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

Generate the Setaria biomass/growth plots

a10.gca.plot = plot_gca('A10')
b100.gca.plot = plot_gca('B100')
r102.gca.plot = plot_gca('R102')
r128.gca.plot = plot_gca('R128')
r133.gca.plot = plot_gca('R133')
r161.gca.plot = plot_gca('R161')
r187.gca.plot = plot_gca('R187')
r20.gca.plot = plot_gca('R20')
r70.gca.plot = plot_gca('R70')
r98.gca.plot = plot_gca('R98')

Supplemental Figure S3

multiplot(a10.gca.plot, b100.gca.plot, r102.gca.plot, r128.gca.plot, 
          r133.gca.plot, r161.gca.plot, r187.gca.plot, r20.gca.plot, 
          r70.gca.plot, r98.gca.plot, cols=3)

agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)

Water use efficiency

In this section we combine data on the volume of water added to each plant per water application (up to the target weight) with fresh-weight biomass data to estimate water use efficiency.

Download water volume data (n = 33,496 water applications).

if (!file.exists('B2_watering_data_phenofront.csv')) {
  download.file('http://files.figshare.com/1845363/B2_watering_data_phenofront.csv',
                'B2_watering_data_phenofront.csv')
}

Read data and format for analysis

water.data = read.table(file="B2_watering_data_phenofront.csv", sep=",", header=TRUE)

Remove the empty pot controls (Barcodes start with zero).

water.data = water.data[grep('000A', water.data$plant.barcode, invert = TRUE),]

Remove snapshots that were not valid waterings (water.amount = -1).

water.data = water.data[water.data$water.amount != -1,]

Add a treatment group column coded in the barcode.

water.data$treatment <- NA
water.data$treatment[grep("AA", water.data$plant.barcode)] <- 100
water.data$treatment[grep("AB", water.data$plant.barcode)] <- 0
water.data$treatment[grep("AC", water.data$plant.barcode)] <- 16
water.data$treatment[grep("AD", water.data$plant.barcode)] <- 33
water.data$treatment[grep("AE", water.data$plant.barcode)] <- 66

Add a plant genotype column coded in the barcode.

water.data$genotype <- NA
water.data$genotype[grep("p1", water.data$plant.barcode)] <- 'A10'
water.data$genotype[grep("p2", water.data$plant.barcode)] <- 'B100'
water.data$genotype[grep("r1", water.data$plant.barcode)] <- 'R20'
water.data$genotype[grep("r2", water.data$plant.barcode)] <- 'R70'
water.data$genotype[grep("r3", water.data$plant.barcode)] <- 'R98'
water.data$genotype[grep("r4", water.data$plant.barcode)] <- 'R102'
water.data$genotype[grep("r5", water.data$plant.barcode)] <- 'R128'
water.data$genotype[grep("r6", water.data$plant.barcode)] <- 'R133'
water.data$genotype[grep("r7", water.data$plant.barcode)] <- 'R161'
water.data$genotype[grep("r8", water.data$plant.barcode)] <- 'R187'

Add a genotype x treatment group column.

water.data$group = paste(water.data$genotype,'-',water.data$treatment,sep='')

Encode the calendar time column as a date-time.

water.data$timestamp = ymd_hms(water.data$timestamp)

Add a column for days after planting since the planting date.

water.data$dap = as.numeric(water.data$timestamp - planting_date)

Some plants were removed from the system but were not inactivated, so remove them.

bad.cars = unique(water.data$car.tag[water.data$water.amount >120 &
                                     water.data$dap > 14])
water.data = water.data[!water.data$car.tag %in% bad.cars,]

Data for Figure 1B

Extract water data for S. viridis for Figure 1B.

sv.water = water.data[(water.data$treatment == 100 | water.data$treatment == 33) &
                       water.data$genotype == 'A10',]

Classify water job as early or later in the day. Early water treatments started at ZT23, later treatments started at ZT8.

sv.water$early.late <- NA  #restrict to the scheduled watering later in the run
sv.water$early.late[hour(sv.water$timestamp) < 7] <- 'ZT23'
sv.water$early.late[hour(sv.water$timestamp) < 15 &
                    hour(sv.water$timestamp) > 11 ] <- 'ZT8'
sv.water = sv.water[!is.na(sv.water$early.late),]
sv.water$dap = day(sv.water$timestamp) + 4
sv.water = sv.water[sv.water$dap > 14,]

Plot the S. viridis water volume data

water.plot = ggplot(sv.water, aes(y = water.amount, x = dap,
                                  group = factor(interaction(early.late,
                                                             treatment,dap)),
                                  fill = factor(interaction(early.late,
                                                            treatment)),
                                  color = factor(interaction(early.late,
                                                             treatment)))) +
                     geom_boxplot(outlier.colour = NULL) +
                     scale_x_continuous("Days after planting") +
                     scale_y_continuous("Water volume (ml)") +
                     theme_bw() +
                     theme(legend.position = c(0.2,0.8),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(fill='Water treatment', color='Water treatment')
water.plot = water.plot + labs(title='Figure 1B')
print(water.plot)

WUE function

Calculate WUE for S. viridis and S. italica

wue.parents = function(water, biomass) {
  # WUE data vectors
  dap.list = c()
  plant.list = c()
  water.list = c()
  biomass.list = c()
  treatment.list = c()
  group.list=c()
  
  # Get unique barcodes for biomass
  barcodes = unique(biomass[(biomass$genotype=='A10' | biomass$genotype=='B100') &
                            (biomass$treatment == 100 |
                             biomass$treatment == 33),]$plant_id)
  for(barcode in barcodes) {
    snapshots = biomass[biomass$plant_id==barcode,]
    snapshots = snapshots[with(snapshots, order(dap)),]
    for(row in 1:nrow(snapshots)) {
      total.water = sum(water[water$dap <= snapshots[row,]$dap &
                              water$plant.barcode == barcode,]$water.amount)
      
      dap.list = c(dap.list, snapshots[row,]$dap)
      plant.list = c(plant.list, barcode)
      water.list = c(water.list, total.water)
      biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
      treatment.list = c(treatment.list, snapshots[row,]$treatment)
      group.list = c(group.list,snapshots[row,]$group)
    }
  }
  
  wue.data = data.frame(plant_id=plant.list,
                        dap=dap.list,
                        water=water.list,
                        biomass=biomass.list,
                        treatment=treatment.list,
                        group=group.list)
  return(wue.data)
}

parents.wue = wue.parents(water.data, vis.data)
parents.wue = parents.wue[parents.wue$water != 0,]

Plot S. viridis and S. italica WUE in mg/mL.

wue.plot = ggplot(parents.wue, aes(x=dap, y=(biomass/water)*1000,
                                   color=factor(group))) +
                  geom_point(size=2.5) +
                  geom_smooth(method="loess", size=1) +
                  scale_x_continuous(name="Days after planting") +
                  scale_y_continuous(lim=c(-2.1, 26),
                                     name="Water-use efficiency (mg/mL)") +
                  theme_bw() +
                  theme(legend.position=c(0.2,0.8),
                        axis.title.x=element_text(face="bold"),
                        axis.title.y=element_text(face="bold")) +
                  labs(color='Genotype-treatment')
wue.plot = wue.plot + labs(title='Figure 4D')
print(wue.plot)

Statistical analysis of WUE

Genotype WUE function

wue = function(water, biomass, genotype) {
  # WUE data vectors
  dap.list = c()
  plant.list = c()
  water.list = c()
  biomass.list = c()
  treatment.list = c()
  
  # Get unique barcodes for biomass
  barcodes = unique(biomass[biomass$genotype==genotype &
                           (biomass$treatment == 100 |
                            biomass$treatment == 33),]$plant_id)
  for(barcode in barcodes) {
    snapshots = biomass[biomass$plant_id==barcode,]
    snapshots = snapshots[with(snapshots, order(dap)),]
    for(row in 1:nrow(snapshots)) {
      total.water = sum(water[water$dap <= snapshots[row,]$dap &
                              water$plant.barcode == barcode,]$water.amount)
      
      dap.list = c(dap.list, snapshots[row,]$dap)
      plant.list = c(plant.list, barcode)
      water.list = c(water.list, total.water)
      biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
      treatment.list = c(treatment.list, snapshots[row,]$treatment)
    }
  }
  
  wue.data = data.frame(plant_id=plant.list,
                        dap=dap.list,
                        water=water.list,
                        biomass=biomass.list,
                        treatment=treatment.list)
  return(wue.data)
}

Analyse pair-wise treatment differences for two-day increments for each genotype.

analyze_wue = function(wue.data) {
  days = c()
  diff.low = c()
  diff.up = c()
  pvals = c()
  for(day in levels(factor(as.integer(wue.data$dap)))) {
    day = as.integer(day)
    control = wue.data[(as.integer(wue.data$dap) == day |
                        as.integer(wue.data$dap) == day + 1) &
                        wue.data$treatment == 100,]
    drought = wue.data[(as.integer(wue.data$dap) == day |
                        as.integer(wue.data$dap) == day + 1) &
                        wue.data$treatment == 33,]
    control.wue = control$biomass / control$water
    drought.wue = drought$biomass / drought$water
    test = t.test(x=control.wue,y=drought.wue)
    days = c(days,day)
    diff.low = c(diff.low, test$conf.int[1])
    diff.up = c(diff.up, test$conf.int[2])
    pvals = c(pvals, test$p.value)
  }
  results = data.frame(dap=as.numeric(days),
                       conf.int.low=diff.low,
                       conf.int.up=diff.up,
                       pvalue=pvals)
  return(results)
}

Analyze WUE for each genotype and aggregate results

wue.results = data.frame()
for(genotype in c('A10', 'B100')) {
  genotype.wue = wue(water.data, vis.data, genotype)
  genotype.wue = genotype.wue[genotype.wue$water != 0,]
  genotype.wue.results = analyze_wue(genotype.wue)
  genotype.wue.results$genotype = genotype
  wue.results = rbind(wue.results, genotype.wue.results)
}

Control for multiple testing by controlling the FDR

qvalues.wue = qvalue(wue.results$pvalue)
wue.results$qvalue = qvalues.wue$qvalues
write.table(wue.results, file='wue.stats.csv', sep=',', row.names=FALSE)
print(wue.results)
##    dap  conf.int.low   conf.int.up       pvalue genotype     qvalue
## 1   11  0.0001290017  6.743361e-04 0.0055453872      A10 0.04066617
## 2   12 -0.0003304552  6.531982e-04 0.5007434856      A10 0.90056523
## 3   13 -0.0004887093  6.645100e-04 0.7558400553      A10 0.95533507
## 4   14 -0.0005521785  8.912840e-04 0.6285589728      A10 0.92188649
## 5   15 -0.0006900247  1.272981e-03 0.5413573872      A10 0.91614327
## 6   16 -0.0013750018  1.352784e-03 0.9863361510      A10 0.98633615
## 7   17 -0.0032500144  1.960050e-03 0.5946773334      A10 0.92188649
## 8   18 -0.0026416988  1.024200e-03 0.3665557006      A10 0.80642254
## 9   19 -0.0025138146  1.198806e-03 0.4687644330      A10 0.90056523
## 10  20 -0.0021111216  2.216305e-03 0.9594528855      A10 0.98633615
## 11  21 -0.0013777865  3.004615e-03 0.4446854076      A10 0.90056523
## 12  22 -0.0012554051  3.263703e-03 0.3590888551      A10 0.80642254
## 13  23 -0.0029737246  4.659146e-03 0.6177633347      A10 0.92188649
## 14  25 -0.0016604803  2.491094e-03 0.6763471555      A10 0.95533507
## 15  26 -0.0016183794  2.333780e-03 0.7068814380      A10 0.95533507
## 16  27 -0.0014579106  1.859317e-03 0.8033499452      A10 0.95533507
## 17  28 -0.0015730725  1.687446e-03 0.9422861297      A10 0.98633615
## 18  29 -0.0016425521  1.397919e-03 0.8687331415      A10 0.98633615
## 19  30 -0.0017235963  1.223038e-03 0.7270761809      A10 0.95533507
## 20  31 -0.0017633197  1.027390e-03 0.5898632801      A10 0.92188649
## 21  32 -0.0018009970  9.249981e-04 0.5116847871      A10 0.90056523
## 22  33 -0.0019234072  2.482920e-03 0.7841395855      A10 0.95533507
## 23  11 -0.0004663250  2.293599e-04 0.4748570604     B100 0.90056523
## 24  12 -0.0003827145  4.190170e-04 0.9249400490     B100 0.98633615
## 25  13 -0.0007208859  1.779908e-04 0.2235222936     B100 0.64019648
## 26  14 -0.0007089761  6.505828e-04 0.9282507680     B100 0.98633615
## 27  15 -0.0013128372  2.973716e-04 0.1995779600     B100 0.62724502
## 28  16 -0.0014265056  1.455777e-03 0.9826344845     B100 0.98633615
## 29  17 -0.0054259745 -3.416845e-04 0.0302482577     B100 0.12099303
## 30  18 -0.0039584923  7.814527e-05 0.0582767730     B100 0.21368150
## 31  19 -0.0053958070 -8.822144e-04 0.0097441669     B100 0.06124905
## 32  20 -0.0049249253 -5.837513e-04 0.0161879273     B100 0.07914098
## 33  21 -0.0062171326 -1.283934e-03 0.0055344094     B100 0.04066617
## 34  22 -0.0056491252 -1.224773e-03 0.0043669943     B100 0.04066617
## 35  23 -0.0061377032 -1.581340e-03 0.0036480639     B100 0.04066617
## 36  25 -0.0056519496 -1.646284e-03 0.0009985777     B100 0.03550789
## 37  26 -0.0050694636 -1.330337e-03 0.0016139948     B100 0.03550789
## 38  27 -0.0041897428 -5.455619e-04 0.0130667633     B100 0.07186720
## 39  28 -0.0036907493 -2.531025e-04 0.0261291938     B100 0.11496845
## 40  29 -0.0033005146  1.619285e-04 0.0734489761     B100 0.24859653
## 41  30 -0.0026887815  6.877754e-04 0.2327987196     B100 0.64019648
## 42  31 -0.0025257356  7.093671e-04 0.2546041527     B100 0.65897545
## 43  32 -0.0024342851  8.292998e-04 0.3164971519     B100 0.77365970
## 44  33 -0.0028628094  2.203841e-03 0.7831321502     B100 0.95533507

Model tiller number

In this section we use linear modeling to estimate tiller number.

Calculate height-width ratio.

vis.data$height_width_ratio = vis.data$height_above_bound / vis.data$extent_x

Download data for manually measured tiller counts 195 random images.

if (!file.exists('200tiller_counts.txt')) {
  download.file('http://files.figshare.com/1845359/200tiller_counts.txt',
                '200tiller_counts.txt')
}

Read manual tiller count data.

tiller.counts = read.table(file="200tiller_counts.txt",sep="\t",header=TRUE)

Format date.

tiller.counts$date = as.POSIXct(paste(paste(tiller.counts$year,
                                            tiller.counts$month,
                                            tiller.counts$day, sep='-'),
                                      paste(tiller.counts$hours,
                                            tiller.counts$minutes,
                                            tiller.counts$seconds, sep=':'), sep=' '))

Merge with VIS snapshot table.

tiller.counts = merge(x=tiller.counts, y=vis.data, by = 'date')

Model tiller counts.

tiller.model.full = lm(ave_tillers ~ fw_biomass + height_above_bound + solidity +
                         extent_x + height_width_ratio, data=tiller.counts)

Variance inflation factor.

vif(tiller.model.full)
##         fw_biomass height_above_bound           solidity 
##          17.338117          22.493335           1.651753 
##           extent_x height_width_ratio 
##          14.767073           3.280364

Drop height and recalculate VIF

tiller.model1 = lm(ave_tillers ~ fw_biomass + solidity + extent_x +
                     height_width_ratio, data=tiller.counts)
vif(tiller.model1)
##         fw_biomass           solidity           extent_x 
##          10.145172           1.393258          11.362312 
## height_width_ratio 
##           1.855525

Drop extent x.

tiller.model2 = lm(ave_tillers ~ fw_biomass + solidity + height_width_ratio,
                   data=tiller.counts)
vif(tiller.model2)
##         fw_biomass           solidity height_width_ratio 
##           1.027223           1.044973           1.039641
summary(tiller.model2)
## 
## Call:
## lm(formula = ave_tillers ~ fw_biomass + solidity + height_width_ratio, 
##     data = tiller.counts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7006 -1.2510 -0.0308  0.8410  5.6311 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5.17731    0.54049   9.579  < 2e-16 ***
## fw_biomass          0.22888    0.01275  17.957  < 2e-16 ***
## solidity           -0.07202    1.87918  -0.038    0.969    
## height_width_ratio -2.16290    0.32335  -6.689 2.42e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.699 on 191 degrees of freedom
## Multiple R-squared:  0.6501, Adjusted R-squared:  0.6446 
## F-statistic: 118.3 on 3 and 191 DF,  p-value: < 2.2e-16

Drop solidity

tiller.model3 = lm(ave_tillers ~ fw_biomass + height_width_ratio, data=tiller.counts)
summary(tiller.model3)
## 
## Call:
## lm(formula = ave_tillers ~ fw_biomass + height_width_ratio, data = tiller.counts)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7003 -1.2525 -0.0264  0.8392  5.6326 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          5.1640     0.4136  12.487  < 2e-16 ***
## fw_biomass           0.2289     0.0126  18.173  < 2e-16 ***
## height_width_ratio  -2.1650     0.3177  -6.815 1.18e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.694 on 192 degrees of freedom
## Multiple R-squared:  0.6501, Adjusted R-squared:  0.6465 
## F-statistic: 178.4 on 2 and 192 DF,  p-value: < 2.2e-16

Plot partial regressions.

avPlots(model = tiller.model3, pch=16, layout=c(2,1), main="Figure 5B")

Download data for manually measured tiller counts for 646 random images.

if (!file.exists('660tiller_counts.txt')) {
  download.file('http://files.figshare.com/1845358/660tiller_counts.txt',
                '660tiller_counts.txt')
}

Predict tiller count for a second set of ~660 randomly selected plants.

tiller660 = read.table(file="660tiller_counts.txt", sep="\t", header=TRUE)
tiller660 = merge(x=tiller660, y=vis.data, by='datetime')
tiller660$predicted_tiller_count = predict.lm(object = tiller.model3,
                                              newdata=tiller660)

Plot manual versus predicted tiller counts.

tiller.plot = ggplot(data=tiller660, aes(x=tiller_count, y=predicted_tiller_count,
                                         color=factor(genotype))) +
                     geom_point(size=2.5) +
                     geom_abline(intercept=0, slope=1) +
                     scale_x_continuous(lim=c(-1,14), "Manual tiller count") +
                     scale_y_continuous(lim=c(-1,14), "Predicted tiller count") +
                     theme_bw() +
                     theme(legend.position=c(0.45,0.95),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold"),
                           legend.direction="horizontal") +
                     labs(color="Genotype")
tiller.plot = tiller.plot + labs(title='Figure 5C')
print(tiller.plot)

Predict tiller counts for the whole data set.

vis.data$tiller_count = predict.lm(object = tiller.model3, newdata = vis.data)

Plot tiller counts for S. viridis and S. italica water treatments 100% and 33% full-capacity.

tc.par.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
                                vis.data$genotype == 'B100') &
                               (vis.data$treatment == 100 |
                                vis.data$treatment == 33),],
                     aes(x=dap, y=tiller_count, color=factor(group))) +
                     geom_point(size=2.5) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(0,12),
                                        name="Estimated tiller count") +
                     theme_bw() +
                     theme(legend.position=c(0.25,0.75),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Genotype-treatment")
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
tc.par.plot = tc.par.plot + labs(title='Figure 5D')
print(tc.par.plot)
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).

Plot estimated tiller counts for all genotypes.

tiller.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
                                    vis.data$treatment == 33,],
                           aes(x=dap, y=tiller_count, color=factor(treatment))) +
                     geom_point(size=1) +
                     geom_smooth(method="loess",size=1) +
                     scale_x_continuous(name="Days after planting") +
                     scale_y_continuous(lim=c(0,12),
                                        name="Estimated tiller count") +
                     theme_bw() +
                     theme(legend.position=c(0.8,0.1),
                           axis.title.x=element_text(face="bold"),
                           axis.title.y=element_text(face="bold")) +
                     labs(color="Treatment") +
                     facet_wrap(~ genotype, ncol = 3)
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 3 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
tiller.facet.plot = tiller.facet.plot + labs(title='Supplemental Figure S4')
print(tiller.facet.plot)
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 3 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).

Trait table

Output a table of all of the traits input and added here. Export all of the VIS traits (except color).

h.table = vis.data
h.table$plant_id = as.character(h.table$plant_id)
h.table$wue = NA
for(r in 1:nrow(h.table)) {
  row = h.table[r,]
  total.water = sum(water.data[water.data$dap <= row$dap & water.data$plant.barcode == row$plant_id,]$water.amount)
  h.table[r,]$wue = row$sv_area / total.water
}
write.table(h.table,file="vis.traits.csv",quote=FALSE,sep=",",row.names=FALSE)